#include "mcmcImageFilters.h"
#include "mcEdgeOp.h"

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////
MCMC_FilterBank::MCMC_FilterBank()
{
	p_filters = NULL;
	p_filtering_result = NULL;
	p_filters_flipped = NULL;

	p_histogram_obs = NULL;
	p_min_value = NULL;
	p_max_value = NULL;
	p_deltas = NULL;
}

MCMC_FilterBank::~MCMC_FilterBank()
{
	if (p_filters != NULL)
		delete []p_filters;
	p_filters = NULL;

	if (p_filters_flipped != NULL)
		delete []p_filters_flipped;
	p_filters_flipped = NULL;

	if (p_filtering_result != NULL)
		delete []p_filtering_result;
	p_filtering_result = NULL;

	FreeHistograms();
}

void MCMC_FilterBank::FreeHistograms(void)
{
	if (p_histogram_obs != NULL)
		delete []p_histogram_obs;
	p_histogram_obs = NULL;

	if (p_min_value != NULL)
		delete []p_min_value;
	p_min_value = NULL;

	if (p_max_value != NULL)
		delete []p_min_value;
	p_min_value = NULL;

	if (p_deltas != NULL)
		delete []p_deltas;
	p_deltas = NULL;
}

void MCMC_FilterBank::Set(Raster<float> &ras_org, const int num)
{
	// Siyuan
	//ras_image = ras_org;
	ras_image.SetDimension(ras_org.nrow, ras_org.ncol);
	ras_org.data.convertTo(ras_image.data, CV_64F);

	num_filter = num;

	if (p_filters != NULL)
		delete []p_filters;
	p_filters = new SparseMatrix<double>[num_filter];

	if (p_filters_flipped != NULL)
		delete []p_filters_flipped;
	p_filters_flipped = new SparseMatrix<double>[num_filter];

	if (p_filtering_result != NULL)
		delete []p_filtering_result;
	p_filtering_result = new Raster<double>[num_filter];

	// pre-defined domain values
	vc_min.SetRowVector(num_filter);
	vc_max.SetRowVector(num_filter);
	vc_predefined_minmax.SetDimension(num_filter,1);
	for (int i=0; i<vc_predefined_minmax.size(); i++)
		vc_predefined_minmax(i) = false;
}

void MCMC_FilterBank::SetOneBank(Raster<double> &ras_filter, const int k, double *pmin, double *pmax)
{
	if (k>=0 && k<NumFilter() && p_filters!=NULL)
	{
		p_filters[k].Transfer(ras_filter);
		p_filters_flipped[k].Transfer(ras_filter);
		p_filters_flipped[k].Flip();

		if (pmin!=NULL && pmax!=NULL)
		{
			vc_min(k) = *pmin;
			vc_max(k) = *pmax;
			vc_predefined_minmax(k) = true;
		}
	}
}

void MCMC_FilterBank::DoFiltering(void)
{
	for (int i=0; i<NumFilter(); i++)
		p_filtering_result[i] = p_filters[i].SparseConv(ras_image);
}


void MCMC_FilterBank::ComputeHistograms(const int cut_off, const bool bcolor)
{
	FreeHistograms();

	p_histogram_obs = new Histogram[NumFilter()];
	p_min_value = new double[NumFilter()];
	p_max_value = new double[NumFilter()];
	p_deltas = new double[NumFilter()];

	int k,i,j;
	for (k=0; k<NumFilter(); k++)
	{
		if (bcolor)
		{
			p_min_value[k] = p_filtering_result[k].GetMin();
			p_max_value[k] = p_filtering_result[k].GetMax();
			if (k==0)
				p_deltas[k] = floor(10.0*max(5.0,(p_max_value[k]-p_min_value[k])/15.0))/10.0;
			else
				p_deltas[k] = floor(10.0*max(1.0,(p_max_value[k]-p_min_value[k])/15.0))/10.0;
		}
		else
		{
			p_min_value[k] = p_filtering_result[k].GetMin();
			p_max_value[k] = p_filtering_result[k].GetMax();

			if (vc_predefined_minmax(k))
			{
				p_min_value[k] = vc_min(k);
				p_max_value[k] = vc_max(k);
			}
			if (k==0)
				p_deltas[k]		= INTENSITY_DELTA;
			else
				p_deltas[k] = floor(10.0*max(1.0,(p_max_value[k]-p_min_value[k])/30.0))/10.0;
		}
		p_histogram_obs[k].Set(p_min_value[k],p_max_value[k],p_deltas[k]);

		for (j=cut_off; j<p_filtering_result[k].rows()-cut_off; j++)
			for (i=cut_off; i<p_filtering_result[k].cols()-cut_off; i++)
				p_histogram_obs[k].AddBin(p_filtering_result[k](j,i),1.0);
		p_histogram_obs[k].Normalize();
	}
}

void MCMC_FilterBank::FreeFilteringResult(void)
{
	if (p_filtering_result != NULL)
		delete []p_filtering_result;
	p_filtering_result = NULL;
}

void MCMC_FilterBank::ComputeHistogram4Mask(mcmcMatrix<int> &mx_mask, const int kcluster,const int cut_off,const bool bcolor)
{
	int k,i,j;
	for (k=0; k<NumFilter(); k++)
	{
		p_min_value[k] = 1e10;
		p_max_value[k] = -1e10;
		// recompute the minimal and the maximal values
		for (j=cut_off; j<p_filtering_result[k].rows()-cut_off; j++)
			for (i=cut_off; i<p_filtering_result[k].cols()-cut_off; i++)
				if (mx_mask(j,i)==kcluster)
				{
					if (p_filtering_result[k](j,i)<p_min_value[k])
						p_min_value[k] = p_filtering_result[k](j,i);
					if (p_filtering_result[k](j,i)>p_max_value[k])
						p_max_value[k] = p_filtering_result[k](j,i);
				}
		if (bcolor)
		{
			if (k==0)
				p_deltas[k] = floor(10.0*max(5.0,(p_max_value[k]-p_min_value[k])/15.0))/10.0;
			else
				p_deltas[k] = floor(10.0*max(1.0,(p_max_value[k]-p_min_value[k])/15.0))/10.0;
		}
		else
		{
			if (vc_predefined_minmax(k))
			{
				p_min_value[k] = vc_min(k);
				p_max_value[k] = vc_max(k);
			}
			if (k==0)
				p_deltas[k]		= INTENSITY_DELTA;
			else
				p_deltas[k] = floor(10.0*max(1.0,(p_max_value[k]-p_min_value[k])/30.0))/10.0;
		}

		p_histogram_obs[k].Set(p_min_value[k],p_max_value[k],p_deltas[k]);
		for (j=cut_off; j<p_filtering_result[k].rows()-cut_off; j++)
			for (i=cut_off; i<p_filtering_result[k].cols()-cut_off; i++)
				if (mx_mask(j,i)==kcluster)
					p_histogram_obs[k].AddBin(p_filtering_result[k](j,i),1.0);
		p_histogram_obs[k].Normalize();
	}
}

void MCMC_FilterBank::ExpandImage(const int row_expansion, const int col_expansion)
{
	Raster<double>	ras_temp;
	int				i,j;

	ras_temp.SetDimension(ras_image.rows()+2*row_expansion, ras_image.cols()+2*col_expansion);
	ras_temp.InitValue(0);

	for (j=row_expansion; j<row_expansion+ras_image.rows(); j++)
		for (i=col_expansion; i<col_expansion+ras_image.cols(); i++)
			ras_temp(j,i) = ras_image(j-row_expansion,i-col_expansion);
	// top
	for (j=0; j<row_expansion; j++)
		for (i=col_expansion; i<col_expansion+ras_image.cols(); i++)
			ras_temp(j,i) = ras_image(row_expansion-1-j,i-col_expansion);
	// bottom
	for (j=row_expansion+ras_image.rows(); j<ras_temp.rows(); j++)
		for (i=col_expansion; i<col_expansion+ras_image.cols(); i++)
			ras_temp(j,i) = ras_image(ras_image.rows()-1-(j-row_expansion-ras_image.rows()),i-col_expansion);
	// left
	for (j=row_expansion; j<row_expansion+ras_image.rows(); j++)
		for (i=0; i<col_expansion; i++)
			ras_temp(j,i) = ras_image(j-row_expansion,col_expansion-1-i);
	// right
	for (j=row_expansion; j<row_expansion+ras_image.rows(); j++)
		for (i=col_expansion+ras_image.cols(); i<ras_temp.cols(); i++)
			ras_temp(j,i) = ras_image(j-row_expansion,ras_image.cols()-1-(i-col_expansion-ras_image.cols()));
	// four corners
	for (j=0; j<row_expansion; j++)
		for (i=0; i<col_expansion; i++)
		{
			ras_temp(j,i) = ras_image(row_expansion-1-j, col_expansion-1-i);
			ras_temp(j,i+col_expansion+ras_image.cols()) = ras_image(row_expansion-1-j, ras_image.cols()-1-i);
			ras_temp(j+row_expansion+ras_image.rows(),i) = ras_image(ras_image.rows()-1-j, col_expansion-1-i);
			ras_temp(j+row_expansion+ras_image.rows(),i+col_expansion+ras_image.cols())
				= ras_image(ras_image.rows()-1-j,ras_image.cols()-1-i);
		}
	ras_image = ras_temp;
}

void MCMC_FilterBank::ShrinkImageBack(const int row_expansion, const int col_expansion)
{
	Raster<double>	ras_temp;
	int				i,j;

	ras_temp.SetDimension(ras_image.rows()-2*row_expansion, ras_image.cols()-2*col_expansion);
	for (j=0; j<ras_image.rows()-2*row_expansion; j++)
		for (i=0; i<ras_image.cols()-2*col_expansion; i++)
			ras_temp(j,i) = ras_image(j+row_expansion,i+col_expansion);

	ras_image = ras_temp;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// filter bank response

// constructions & destructions
MCMC_FilterResp::MCMC_FilterResp()
{
	vc_nbc.SetRowVector(4);
	vc_nbc(0) = -1;   vc_nbc(1) = 0;		vc_nbc(2) = 0;		vc_nbc(3) = 1;
	vc_nbr.SetRowVector(4);
	vc_nbr(0) = 0;    vc_nbr(1) = -1;		vc_nbr(2) = 1;		vc_nbr(3) = 0;
}

MCMC_FilterResp::~MCMC_FilterResp()
{
}

void MCMC_FilterResp::Set(Raster_Color<uchar> &raster1, const bool bcolor1)
{
	int		i,j;

	bcolor = bcolor1;

	// color image
	image_original = raster1;
	image_original.CreateLUVspace();

	// intensity
	ras_intensity.SetDimension(image_original.rows(),image_original.cols());
	for (j=0; j<ras_intensity.rows(); j++)
		for (i=0; i<ras_intensity.cols(); i++)
			ras_intensity(j,i) = image_original.gray(j,i);

	min_L = 1e10; max_L = -1e10;
	min_U = 1e10; max_U = -1e10;
	min_V = 1e10; max_V = -1e10;

	// copy L,U,V out
	ras_L.SetDimension(ras_intensity.rows(),ras_intensity.cols());
	ras_U.SetDimension(ras_intensity.rows(),ras_intensity.cols());
	ras_V.SetDimension(ras_intensity.rows(),ras_intensity.cols());
	for (j=0; j<ras_intensity.rows(); j++)
		for (i=0; i<ras_intensity.cols(); i++)
		{
			ras_L(j,i) = image_original.L(j,i);
			if (ras_L(j,i)<min_L)	min_L = ras_L(j,i);
			if (ras_L(j,i)>max_L)	max_L = ras_L(j,i);

			ras_U(j,i) = image_original.U(j,i);
			if (ras_U(j,i)<min_U)	min_U = ras_U(j,i);
			if (ras_U(j,i)>max_U)	max_U = ras_U(j,i);

			ras_V(j,i) = image_original.V(j,i);
			if (ras_V(j,i)<min_V)	min_V = ras_V(j,i);
			if (ras_V(j,i)>max_V)	max_V = ras_V(j,i);
		}

	// get edge pixles on this image
	Raster<double>	Ix,Iy;
	double dfactor;
	
	if (bcolor)
	{
		MCMC_EdgeOp_Color  edge_color;

		edge_color = image_original;
		// Siyuan
		//ras_edge_g = edge_color.Canny(1.5, 0.05, &Ix, &Iy);
		Raster<double> ras_edge_g_double;
		ras_edge_g_double = edge_color.Canny(1.5, 0.05, &Ix, &Iy);
		ras_edge_g.SetDimension(ras_edge_g_double.nrow, ras_edge_g_double.ncol);
		ras_edge_g_double.data.convertTo(ras_edge_g.data, CV_32F);

		dfactor = 6.0;
	}
	else
	{
		MCMC_EdgeOp edge_grey;
		
		// Siyuan
		//edge_grey = ras_intensity;
		//ras_edge_g = edge_grey.Canny(2.0, 0.05, &Ix, &Iy);
		edge_grey.SetDimension(ras_intensity.nrow, ras_intensity.ncol);
		ras_intensity.data.convertTo(edge_grey.data, CV_64F);

		Raster<double> ras_edge_g_double;
		ras_edge_g_double = edge_grey.Canny(2.0, 0.05, &Ix, &Iy);
		ras_edge_g.SetDimension(ras_edge_g_double.nrow, ras_edge_g_double.ncol);
		ras_edge_g_double.data.convertTo(ras_edge_g.data, CV_32F);

		dfactor = 3.0;
	}

	ras_edge_theta.SetDimension(ras_edge_g.rows(), ras_edge_g.cols());
	for (j=0; j<ras_edge_g.rows(); j++)
		for (i=0; i<ras_edge_g.cols(); i++)
		{
			ras_edge_g(j,i) = min(ras_edge_g(j,i)/dfactor, 5.0);
			ras_edge_theta(j,i) = atan2(Iy(j,i), Ix(j,i));
		}

	// get the DC components and High frequency components of the image
	// use a Gaussian kernel to blur the image to subtract DC components
	mcmcMatrix<float> mx_gauss_1d,mx_dest;

	mx_gauss_1d.SetGaussian(1,3.0);
	mx_dest = ras_intensity.ConvHor(mx_gauss_1d);
	mx_dest = mx_dest.ConvVer(mx_gauss_1d.T());
	for (int k=0; k<2; k++)
	{
		mx_dest = mx_dest.ConvHor(mx_gauss_1d);
		mx_dest = mx_dest.ConvVer(mx_gauss_1d.T());
	}
	ras_DC  = mx_dest;

	// subtract the blurred image from the original image
	ras_HF	= ras_intensity-ras_DC;
	// initilize the region layer, residual image
	ras_region_layer = ras_intensity;
	ras_region_rec.SetDimension(rows(),cols());
	ras_region_rec.InitValue(0);
	ras_region_residue.SetDimension(rows(),cols());
	ras_region_residue.InitValue(0);
	ras_image_reconstruction.SetDimension(rows(),cols());
	ras_image_reconstruction.InitValue(0);

	// integral image
	ComputeIntegralImage();
}

void MCMC_FilterResp::ComputeIntegralImage(void)
{
	Raster<uchar> ras_temp;

	CopyRaster(ras_temp, ras_intensity);

	itg_image.Prepare(ras_temp);
}


void MCMC_FilterResp::ComputeHistogramAsIid(RF3D &rf_iid,mcmcMatrix<int> &mx_mask, const int kcluster)
{
	int i,j;

	for (j=0; j<image_original.rows(); j++)
		for (i=0; i<image_original.cols(); i++)
			if (mx_mask(j,i)==kcluster)
			{
				rf_iid.AddOneBin(image_original.L(j,i),image_original.U(j,i),image_original.V(j,i),1.0);
			}
}

void MCMC_FilterResp::ComputeHistogramAsIid(RF &rf_iid,mcmcMatrix<int> &mx_mask, const int kcluster)
{
	int i,j;

	for (j=0; j<image_original.rows(); j++)
		for (i=0; i<image_original.cols(); i++)
			if (mx_mask(j,i)==kcluster)
			{
				rf_iid.AddOneBin(ras_intensity(j,i),1.0);
			}
}


int MCMC_FilterResp::TotoalNumFilterResult(void)
{
	if (bcolor)
		return		(filter_bank_L.NumFilter()+	filter_bank_U.NumFilter()+
					filter_bank_V.NumFilter());
	else
		return		filter_bank_intensity.NumFilter();
}

// get filtering result
Raster<double>* MCMC_FilterResp::GetFirstResult(double &delta)
{
	which_result = 0;

	if (bcolor)
	{
		delta = filter_bank_L.p_deltas[0];
		return &(filter_bank_L.p_filtering_result[0]);
	}
	else
	{
		delta = filter_bank_intensity.p_deltas[0];
		return &(filter_bank_intensity.p_filtering_result[0]);
	}
}

// get filtering result
Raster<double>* MCMC_FilterResp::GetNextResult(double &delta)
{
	which_result++;

	int ntemp=which_result;

	if (ntemp<0)
		return NULL;

	if (bcolor)
	{
		if (ntemp<filter_bank_L.NumFilter())
		{
			delta = filter_bank_L.p_deltas[ntemp];
			return &(filter_bank_L.p_filtering_result[ntemp]);
		}
		ntemp = ntemp-filter_bank_L.NumFilter();

		if (ntemp<filter_bank_U.NumFilter())
		{
			delta = filter_bank_U.p_deltas[ntemp];
			return &(filter_bank_U.p_filtering_result[ntemp]);
		}
		ntemp = ntemp-filter_bank_U.NumFilter();

		if (ntemp<filter_bank_V.NumFilter())
		{
			delta = filter_bank_V.p_deltas[ntemp];
			return &(filter_bank_V.p_filtering_result[ntemp]);
		}
		ntemp = ntemp-filter_bank_V.NumFilter();
	}
	else
	{
		if (ntemp<filter_bank_intensity.NumFilter())
		{
			delta = filter_bank_intensity.p_deltas[ntemp];
			return &(filter_bank_intensity.p_filtering_result[ntemp]);
		}
	}

	return NULL;
}

void MCMC_FilterResp::SetAllFilterBanksForClustering()
{
	Raster<double>	ras_filters[16],ras_LoG[3],ras_sigma;
	int				i,num_filter=2;

	ras_sigma.SetDimension(1,1);
	ras_sigma(0) = 1.0f;

	// filters
	ras_filters[0].SetDimension(3,3);
	ras_filters[0].InitValue(0.0f);
	ras_filters[0](0,0) = 1.0f; ras_filters[0](0,1) = 1.0f; ras_filters[0](0,2) = 1.0f;
	ras_filters[0](2,0) = -1.0f; ras_filters[0](2,1) = -1.0f; ras_filters[0](2,2) = -1.0f;

	ras_filters[1].SetDimension(3,3);
	ras_filters[1].InitValue(0.0f);
	ras_filters[1](0,0) = -1.0f; ras_filters[1](0,2) = 1.0f;
	ras_filters[1](1,0) = -1.0f; ras_filters[1](1,2) = 1.0f;
	ras_filters[1](2,0) = -1.0f; ras_filters[1](2,2) = 1.0f;

	// Lapaliacian of Gaussian
	ras_LoG[0].SetDimension(3,3);
	ras_LoG[0].InitValue(0.0f);
	ras_LoG[0](0,1) = 1.0f;	ras_LoG[0](1,0) = 1.0f;
	ras_LoG[0](0,2) = 1.0f;	ras_LoG[0](2,0) = 1.0f;
	ras_LoG[0](1,1) = -4.0f;

	ras_LoG[1].SetDimension(5,5);
	ras_LoG[1].InitValue(0.0f);
	ras_LoG[1](0,2)=-1;
	ras_LoG[1](1,1)=-1; ras_LoG[1](1,2)=-2; ras_LoG[1](1,3)=-1;
	ras_LoG[1](2,0)=-1; ras_LoG[1](2,1)=-2; ras_LoG[1](2,2)=16; ras_LoG[1](2,3)=-2; ras_LoG[1](2,4)=-1;
	ras_LoG[1](3,1)=-1; ras_LoG[1](3,2)=-2; ras_LoG[1](3,3)=-1;
	ras_LoG[1](4,2)=-1;

	ras_LoG[2].SetDimension(5,5);
	ras_LoG[2].InitValue(0.0f);
	ras_LoG[2](0,2) = 1.0f;
	ras_LoG[2](2,0) = 1.0f; ras_LoG[2](2,2) = -4.0f; ras_LoG[2](2,4) = 1.0f;
	ras_LoG[2](4,2) = 1.0f;

	cut_off = 0;

	if (bcolor)
	{
		// in L
		filter_bank_L.Set(ras_L,4+num_filter);
		filter_bank_L.SetOneBank(ras_sigma,0);
		filter_bank_L.SetOneBank(ras_LoG[0],1);
		filter_bank_L.SetOneBank(ras_LoG[1],2);
		filter_bank_L.SetOneBank(ras_LoG[2],3);

		for (i=0; i<num_filter; i++)
			filter_bank_L.SetOneBank(ras_filters[i],i+4);

		filter_bank_L.DoFiltering();
		filter_bank_L.ComputeHistograms(cut_off);

		// in U
		filter_bank_U.Set(ras_U,3);
		filter_bank_U.SetOneBank(ras_sigma,0);
		filter_bank_U.SetOneBank(ras_LoG[0],1);
		filter_bank_U.SetOneBank(ras_LoG[2],2);


		filter_bank_U.DoFiltering();
		filter_bank_U.ComputeHistograms(cut_off);

		// in V
		filter_bank_V.Set(ras_V,3);
		filter_bank_V.SetOneBank(ras_sigma,0);
		filter_bank_V.SetOneBank(ras_LoG[0],1);
		filter_bank_V.SetOneBank(ras_LoG[2],2);

		filter_bank_V.DoFiltering();
		filter_bank_V.ComputeHistograms(cut_off);
	}
	else
	{
		filter_bank_intensity.Set(ras_intensity,4+num_filter);
		filter_bank_intensity.SetOneBank(ras_sigma,0);
		filter_bank_intensity.SetOneBank(ras_LoG[0],1);
		filter_bank_intensity.SetOneBank(ras_LoG[1],2);
		filter_bank_intensity.SetOneBank(ras_LoG[2],3);

		for (i=0; i<num_filter; i++)
			filter_bank_intensity.SetOneBank(ras_filters[i],i+4);

		filter_bank_intensity.DoFiltering();
		filter_bank_intensity.ComputeHistograms(cut_off,bcolor);
	}
}


double GetRasterPositiveSum(const Raster<double> &filter)
{
	double sum=0.0;

	for (int i=0; i<filter.size(); i++)
		if (filter.Data(i)>0)
			sum += filter.Data(i);
	return sum;
}

double GetRasterNegtiveSum(const Raster<double> &filter)
{
	double sum=0.0;

	for (int i=0; i<filter.size(); i++)
		if (filter.Data(i)<0)
			sum += filter.Data(i);
	return sum;
}


void MCMC_FilterResp::SetAllFilterBanks()
{
	Raster<double>	ras_gradient[2],ras_LoG[3],ras_sigma;
	int				i,k,j,num_LoG=1,num_gabor=4,num_grad=2;
	mcmcVector<double>	vc_sigmas(1),vc_thetas(4);
	Raster<double>	ras_gabors[16];
	int				max_size_filter;

	vc_sigmas(0) = 1.5f; //vc_sigmas(1) = 2.5f;
	vc_thetas(0) = 0.0f; vc_thetas(1) = 90.0f; vc_thetas(2) = 45.0f; vc_thetas(3) = 135.0f;
	ras_sigma.SetDimension(1,1);
	ras_sigma(0) = 1.0f;

	ras_sigma.SetDimension(1,1);
	ras_sigma(0) = 1.0f;

	// Lapaliacian of Gaussian
	ras_LoG[0].SetDimension(3,3);
	ras_LoG[0].InitValue(0.0f);
	ras_LoG[0](0,1) = 0.125;	ras_LoG[0](1,0) = 0.125;
	ras_LoG[0](1,2) = 0.125;	ras_LoG[0](2,1) = 0.125;
	ras_LoG[0](1,1) = -0.5f;

	// filters
	ras_gradient[0].SetDimension(3,3);
	ras_gradient[0].InitValue(0.0f);
	ras_gradient[0](0,0) = 1.0/6.0; ras_gradient[0](0,1) = 1.0/6.0; ras_gradient[0](0,2) = 1.0/6.0;
	ras_gradient[0](2,0) = -1.0/6.0; ras_gradient[0](2,1) = -1.0/6.0; ras_gradient[0](2,2) = -1.0/6.0;

	ras_gradient[1].SetDimension(3,3);
	ras_gradient[1].InitValue(0.0f);
	ras_gradient[1](0,0) = -1.0/6.0; ras_gradient[1](0,2) = 1.0/6.0;
	ras_gradient[1](1,0) = -1.0/6.0; ras_gradient[1](1,2) = 1.0/6.0;
	ras_gradient[1](2,0) = -1.0/6.0; ras_gradient[1](2,2) = 1.0/6.0;

	max_size_filter = 3;
	// Gabor Filters
	k = 0;
	for (j=0; j<vc_sigmas.size(); j++)
		for (i=0; i<vc_thetas.size(); i++)
		{
			ras_gabors[k].SetGabor(1.0,vc_sigmas(j),vc_thetas(i));
			double a=ras_gabors[k].Sum();
			ras_gabors[k] = ras_gabors[k]-a/ras_gabors[k].size();

			// get the size of the biggest filter
			if (ras_gabors[k].rows()>max_size_filter)
				max_size_filter = ras_gabors[k].rows();
			if (ras_gabors[k].cols()>max_size_filter)
				max_size_filter = ras_gabors[k].cols();

			k++;
		}

	cut_off=(int)ceil(((double)max_size_filter)/2.0)+1;

	if (bcolor)
	{
		// color model
		// in L
		filter_bank_L.Set(ras_L,1+num_grad+num_LoG+num_gabor);
		filter_bank_L.SetOneBank(ras_sigma,0);
		for (i=0; i<num_grad; i++)
			filter_bank_L.SetOneBank(ras_gradient[i],i+1);
		for (i=0; i<num_LoG; i++)
			filter_bank_L.SetOneBank(ras_LoG[i],i+1+num_grad);
		for (i=0; i<num_gabor; i++)
			filter_bank_L.SetOneBank(ras_gabors[i],i+1+num_grad+num_LoG);

		filter_bank_L.DoFiltering();
		filter_bank_L.ComputeHistograms(cut_off);

		// in U
		filter_bank_U.Set(ras_U,1+num_grad+num_LoG);
		filter_bank_U.SetOneBank(ras_sigma,0);
		for (i=0; i<num_grad; i++)
			filter_bank_U.SetOneBank(ras_gradient[i],i+1);
		for (i=0; i<num_LoG; i++)
			filter_bank_U.SetOneBank(ras_LoG[i],i+3);

		filter_bank_U.DoFiltering();
		filter_bank_U.ComputeHistograms(cut_off);

		// in V
		filter_bank_V.Set(ras_V,1);
		filter_bank_V.SetOneBank(ras_sigma,0);

		filter_bank_V.DoFiltering();
		filter_bank_V.ComputeHistograms(cut_off);
	}
	else
	{
		double min_value,max_value;

		filter_bank_intensity.Set(ras_intensity,1+num_grad+num_LoG+num_gabor);
		min_value = INTENSITY_LOW; max_value=INTENSITY_HIGH;
		filter_bank_intensity.SetOneBank(ras_sigma,0,&min_value,&max_value);

		for (i=0; i<num_grad; i++)
		{
			min_value = floor(GetRasterNegtiveSum(ras_gradient[i])*255.0*10.0)/10.0;
			max_value = ceil(GetRasterPositiveSum(ras_gradient[i])*255.0*10.0)/10.0;
			filter_bank_intensity.SetOneBank(ras_gradient[i],i+1,&min_value,&max_value);
		}
		for (i=0; i<num_LoG; i++)
		{
			min_value = floor(GetRasterNegtiveSum(ras_LoG[i])*255.0*10.0)/10.0;
			max_value = ceil(GetRasterPositiveSum(ras_LoG[i])*255.0*10.0)/10.0;
			filter_bank_intensity.SetOneBank(ras_LoG[i],i+1+num_grad,&min_value,&max_value);
		}
		for (i=0; i<num_gabor; i++)
		{
			min_value = floor(GetRasterNegtiveSum(ras_gabors[i])*255.0*10.0)/10.0;
			max_value = ceil(GetRasterPositiveSum(ras_gabors[i])*255.0*10.0)/10.0;
			filter_bank_intensity.SetOneBank(ras_gabors[i],i+1+num_grad+num_LoG,&min_value,&max_value);
		}

		filter_bank_intensity.DoFiltering();
		filter_bank_intensity.ComputeHistograms(cut_off,bcolor);
	}
}

void MCMC_FilterResp::ComputeHistogram(void)
{
	if (bcolor)
	{
		filter_bank_L.ComputeHistograms(cut_off,true);
		filter_bank_U.ComputeHistograms(cut_off,true);
		filter_bank_V.ComputeHistograms(cut_off,true);
	}
	else
		filter_bank_intensity.ComputeHistograms(cut_off,false);
}


void MCMC_FilterResp::ComputeHistogram4Mask(mcmcMatrix<int> &mx_mask, const int kcluster)
{
	if (bcolor)
	{
		filter_bank_L.ComputeHistogram4Mask(mx_mask,kcluster,cut_off);
		filter_bank_U.ComputeHistogram4Mask(mx_mask,kcluster,cut_off);
		filter_bank_V.ComputeHistogram4Mask(mx_mask,kcluster,cut_off);
	}
	else
		filter_bank_intensity.ComputeHistogram4Mask(mx_mask,kcluster,cut_off,bcolor);
}


void MCMC_FilterResp::ExpandImage(const int row_expansion, const int col_expansion)
{
	if (bcolor)
	{
		filter_bank_L.ExpandImage(row_expansion,col_expansion);
		filter_bank_U.ExpandImage(row_expansion,col_expansion);
		filter_bank_V.ExpandImage(row_expansion,col_expansion);
	}
	else
		filter_bank_intensity.ExpandImage(row_expansion,col_expansion);
}


void MCMC_FilterResp::ShrinkImageBack(const int row_expansion, const int col_expansion)
{
	if (bcolor)
	{
		filter_bank_L.ShrinkImageBack(row_expansion,col_expansion);
		filter_bank_U.ShrinkImageBack(row_expansion,col_expansion);
		filter_bank_V.ShrinkImageBack(row_expansion,col_expansion);
	}
	else
		filter_bank_intensity.ShrinkImageBack(row_expansion,col_expansion);
}

void MCMC_FilterResp::DoFiltering(void)
{
	if (bcolor)
	{
		filter_bank_L.DoFiltering();
		filter_bank_U.DoFiltering();
		filter_bank_V.DoFiltering();
	}
	else
		filter_bank_intensity.DoFiltering();
}

void MCMC_FilterResp::FilteringOnNewImage(Raster<uchar> &image_intensity)
{
	ras_intensity.SetDimension(image_intensity.rows(),image_intensity.cols());
	for (int j = 0; j < ras_intensity.rows(); j++){
		for (int i = 0; i < ras_intensity.cols(); i++){
			ras_intensity(j, i) = (double)image_intensity(j, i);
		}
	}

	// Siyuan
	//filter_bank_intensity.ras_image = ras_intensity;
	filter_bank_intensity.ras_image.SetDimension(ras_intensity.nrow, ras_intensity.ncol);
	ras_intensity.data.convertTo(filter_bank_intensity.ras_image.data, CV_64F);

	filter_bank_intensity.DoFiltering();
}
